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Abstract 

<^ ' We present two algorithms for the computation of the Kalman form of a linear control system. 

The first one is based on the technique developed by Keller-Gehrig for the computation of 
the characteristic polynomial. The cost is a logarithmic number of matrix multiplications. To 
our knowledge, this improves the best previously known algebraic complexity by an order of 
magnitude. Then we also present a cubic algorithm proven to be more efficient in practice. 

O ■ 



■ 1 Introduction 

H 

This report is a continuation of a collaboration of the first two authors on the algorithmic simi- 
larities between the computation of the Kalman form and of the characteristic polynomial. This 
collaboration led to [DR05, Theorem 2]. We report here an improvement of this last result based 
on a remark by the third author. 

For a definition of the Kalman form of a linear control system, see [Kal61 , Theorem 1]. 
In this report we show how to adapt the branching algorithm of Keller-Gehrig [KG85, §5] 
(computing the characteristic polynomial) to compute the Kalman form. This implies an algebraic 
time complexity of 0(n"logn). Now, the discussion of [DPW05, §2] shows that a cubic algorithm, 
luk, is more efficient in practice for the computation of the characteristic polynomial. Therefore, 
we adapt it to the computation of the Kalman form. 

The outline of this report is the following : in section 2 we define the compressed Kyrlov 
matrix. It will help to describe the adaptation of Keller-Gehrig's algorithm to the computation of the 
Kalman form. In section 3, we recall Keller-Gehrig's algorithm. Section 4 presents the main result 
of this report, on the time complexity of the computation of the Kalman form. Lastly, we give a 
full description two algorithms to compute the Kalman form. The first one precises the operations 
used in section 4 to achieve the complexity and improves the constant hidden in the 0{) by saving 
operations. The second is based on another algorithm for the characteristic polynomial, that does 
not achieve the same algebraic complexity, but appears to be faster in practice. 
We will denote by u> the exponent in the complexity of the matrix multiplication. 
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2 The compressed Krylov matrix 



Let A and B be two matrices of dimension respectively n x n and n x m. Consider the n x (mri) 
Krylov matrix if generated by the m column vectors of B and their iterates with the matrix A : 

K=[h \ ... | A^h | ... | 6 ro | ... | A n ~ 1 b m } 

Let r be the rank of K. r > rank(B). Let us form the n x r non-singular matrix ~K by picking 
the first r linearly independent columns of K. 

Definition 2.1 . ~K is the compressed Krylov matrix of B relatively to A. 

If a column vector A k b 3 is linearly dependent with the previous column vectors, then any vector 
A l bj, I > k will also be linearly dependent. Consequently the matrix ~K has the form : 

K=[by \ ... | | ... | b m | ... | A dm ~ x b m ] (1) 

for some such that < di < n - 1 and YnLi d i = r - 

The order in the choice of the independent column vectors (from the left to the right) can also 
be interpreted in terms of lexicographical order on the sequence (di). Following Storjohann [?], 
we can therefore also define the compressed Krylov matrix as follows : 

Definition 2.2. The compressed Krylov matrix of B relatively to A is a matrix of the form 

[ h | ... | A^h | ... | b m | ... | A Am ~ x b m ] 

of rank r, such that the sequence (di) is lexicographically maximal. 

The next section will present an algorithm to compute this compressed Krylov matrix. 

3 Keller-Gehrig's algorithm 

The selection of the linearly independent columns, starting from left to right, can be done by 
a gaussian elimination. A block elimination is mandatory to reduce the algebraic complexity to 
matrix multiplication. For this task, Keller-Gehrig first introduced in [KG85, §4] an algorithm called 
"step form elimination". The more recent litterature replaced it by the row echelon elimination (for 
example in [CBS97]). We showed in [DPW05] that the LQUP elimination (defined in [IMH82]) of 
~K T could also be used (algorithm 3.1). This last algorithm simply returns the submatrix formed 
by the first independent column vectors of the input matrix form left to right. 



Algorithm 3.1 ColReducedForm 

Require: A a m x n matrix of rank r (m,n> r) over a field 

Ensure: A' am x r matrix formed by r linearly independent columns of A 

1 : (L, Q, U, P, r) = LQUP(A T ) (r = rank(A)) 
2: return ([/ r 0](Q T ,4 T )) T 



Thus a straightforward algorithm to compute ~K would be to run algorithm 3.1 on the matrix 
K. But cost of the computation of K is prohibitive (n 3 coefficients and 0(n 4 ) arithmetic opera- 
tions with standard matrix product). Hence, the elimination process must be combined within the 
building of the matrix. 

The computation of the iterates can rely on matrix multiplication, by computing the [log 2 (n)] 
following powers of A : 

AA 2 ,...,A 2 \^ no92< " ,W 
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Thus the following scheme, 

Vo = [bj] 

V l+1 = [V l \A T V l ] { ' 

where the matrix V t has 2 i columns, computes every iterates of bj in O(n"logn) operations. 

One elimination is performed after each application of A 2 \ to discard the linearly dependent 
iterates for the next iteration step. Moreover if a vector bj has only k < 2 l linearly independent 
iterates, one can stop the computation of its iterates. Therefore, the scheme (2) will only be 
applied on the block iterates of size 2\ 

From these remarks, we can now present Keller-Gehrig's algorithm. Although is was initially 
designed for the computation of the characteristic polynomial, we prefer to show it in a more 
general setting : the computation of the compressed Krylov matrix. Afterwards, we will show that 
the computation of the characteristic polynomial is a specialization of this algorithm with B = I n 
and that the recover of its coefficients is straightforward. 

Algorithm 3.2 Compressed Krylov Matrix [Keller-Gehrig] 
Require: Aanx n matrix over a field, B, a n x m matrix 
Ensure: (K,r) as in (1) 
1 

2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 



i = 

Vo = B = (y ,i,v ,2,...,v , m ) 

C = A 

while (3k, Vfe has 2* columns) do 
for all j do 

if ( Vij has strictly less than 2* columns ) then 

Wj = 



14 
15 
16 
17 



else 

Wj = [V iij \CV i ,j] 
end if 
end for 

W=[W 1 \...\W n ] 

Vi+i = ColReducedForm(VF) remember r = rank(W) {V l+ i = . . . |^+i,n] where 

Vi+ij are the remaining vectors of Wj in V i+1 } 
C = CxC 
i = i + l 
end while 
return (Vi,r) 



Theorem 3.1 (Keller-Gehrig). Suppose m = 0(n). The compressed Krylov matrix of B relatively 
to A can be computed in 0(n u logn) field operations. 

Proof. Algorithm 3.2 satisfies the statement (cf [KG85]). □ 

Property 3.1 . Let ~K be the compressed Krylov matrix of the identity matix relatively to A. The 
matrix ~K 1 A~K has the Hessenberg poly cyclic form : it is block upper triangular, with companion 
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blocks on its diagonal, and the upper blocks are zero except on their last column. 



K 1 AK = 






* 


1 


* 




* 




1 * 






* 


1 


* 




* 




1 * 



(3) 



Corollary 3.1 (Keller-Gehrig). The characteristic polynomial of A can be computed in O(n w logn) 

field operation. 

Proof. The characteristic polynomial of the shifted Hessenberg form (3) is the product of the 
polynomials associated to the companion blocks on its diagonal. And since determinants are 
invariants under similarity transformations, it equals the characteristic polynomial of A. □ 

4 Computation of the Kalman form 

Theorem 4.1 recalls the definition of the Kalman form of two matrices A and B. 

Theorem 4.1. Let A and B be two matrices of dimension respectively n x n and n x m. Let r 
be the dimension of Span(B, AB, . . . , A n ~ 1 B). There exist a non singular matrix T of dimension 
nxn such that 



T -1 AT 



' H 


X ' 




' B 1 ' 





Y 


1 






= T^B 



where H and B\ are respectively r x r et r x m. 

The main result of this report is the following result, based on an idea by the third author. 

Theorem 4.2. Let V be compressed Krylov matrix of B respectively to A. Complete V into a 
basis T of K n by adding n — r columns at the end of V. Then T satisfies the definition of the 
Kalman form of A and B. 



Proof. The matrix V satisfy the relation 



AV = VH 



[V\W]. 



where H is r x r and has the Hessenberg polycylic form (3). Let us note T 
Now 

AT = [ AV | AW ] 

H X 


Lastly, V is a basis of Span(B, AB, A n B). Therefore each column of B is a linear combi- 
nation of the colmuns of V : 

" B x 



B 







□ 



Corollary 4.1 . The Kalman form of A and B can be computed in 0(n u logn). 
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Proof. Applying theorem 3.1 , there only remains to show how to complete V into T in O(n w logn). 
The idea is to complete V in its triangularized form. One computes the LUP factorization of V T : 



V T = [L}[ U x U 2 }P 



Then replace [UiU 2 ] by 
simply corresponds to set 



Id 



and [L] by 



L 
Id 



to get a n x n non singular matrix. This 



T = 



K 


P T 











In — r 





It only costs 0(n w ) field operations to recover the whole Kalman form (blocks H,X,Y and 
Bi), using for example matrix multiplications and matrix inversions. See section 5.2 for more 
details. □ 

This last result improves the algebraic time complexity for the computation of the Kalman form 
given in [DR05, Theorem 2] by an order of magnitude. 



5 Algorithms into practice 

The goal of the previous section was to establish the time complexity estimate and we therefore 
only sketched the algorithms involved. We will now focus more precisely on the operations so as 
to reduce the consant hiden in the 0{) notation. 



5.1 Improvements on Keller-Gehrig's algorithm 



The first improvement concerns the recover of the Hessenberg polycyclic form 3, once the com- 
pressed Krylov matrix is computed. In [KG85] Keller-Gehrig simply suggests to compute the prod- 
uct K~ X AK. This implies 4.66n 3 additional field operation. We propose here to reduce this cost 
to 4>n 2 , where <fi is the number of blocks in the Hessenberg form. This technique was presented in 
[DPW05]. We recall and extend it here for the recovery of the whole Hessenberg polycylic form. 

First consider the case where the n first iterates of only one vector v are linearly independent. 

Let K = [v\Av\ . . . \A n v\. The last column is the first which is linearly dependent with the 
previous. Let P{X) = X n - Y17=o m i xi represent this dependency (the minimal polynomial of 
this vector relatively to A). Again consider the LUP factorization of K T . Let X n+ i denote the 
n + 1th row of the matrix X and Xi,,, n be the block of the first n rows of X. Then we have 



Therefore 



Kl +1 = (A n v) T = (J2 m^vf = [ m ... m„-i Ki^k 



»=o 



L n+ i = [m ... m n _i ]£i... n . 



And the coefficients m, can be recovered as the solution of a triangular system. 
Now, one easily check that 



K^. n AKi...n 





1 



1 



too 

TOi 



This companion matrix is the Hessenberg polycyclic matrix to be computed. 

In the situation of Keller-Gehrig's algorithm, the linear depencies also involve iterates of other 
vectors. However, the LQUP factorization will play a similar role than the previous LUP and makes 
it possible to recover the whole vector coefficients of the linear dependency. 
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Figure 1 : LQUP factorization of 2 blocks of iterates 




Figure 2: Recover of the coefficients of the first linear dependency 

We show in figure 1 the case of two blocks of iterates. The first linear dependency relation (for 
A dl vi) is done as previously (see figure 2). 

Now for the second block, the first linearly dependent vector A d2 v 2 satisfies a relation of the 
type : 

d.2-1 

A d * V2 = ]T MS + ^ A * Vl 

i=0 i=0 

The vector of coefficients (3 = [fa] and 7 = [7*] can be obtained by solving the following system 
shown in figure 3. 




Figure 3: Recover of the coefficients of the second linear dependency 
There only remains to build the Hessenberg polycyclic matrix from these vectors : 
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H = 





1 



a 
at 



7o 
7i 



7A-1 

A) 
ft 



1 /?d 2 -i 

This technique can be applied to every block of iterates. Therefore the Hessenberg polycyclic 
matrix can be recovered by as many triangular system resolutions as the number of blocks. 



5.2 The main algorithm 

We have seen in section 4 how to compute the matrix T (simply T 





P T 











I n — r 





)■ 



Section 5.1 showed how to compute H. There only remains to show how to compute the matrices 
X, Y and Bi and we will be done. 

We recall (from the proof of theorem 4.2) that 



AT 



AV 


AP T 











I n — r 





= T 



H X 
Y 



Then X and Y satisfy T 



We have 



Now 



therefore 



AP 1 



. Let us write 



A' = PAP T = 



A' A' 

^11 A 12 

A. 



x 2l 



A' 

^22 



'X' 


= P T A' 





= p t 


'A' 

A \2 


Y 




I n — r 




A 1 
.22. 



T = P 1 



~ul " 




L T " 


_ L '2 1 n—r_ 




In — r 



~ul o " 




L T " 




_ ^2 ^n-r_ 




I n — r 





X 
Y 

L T X 
Y 



And the system 



has the following solution 



U{ L 1 X 
UTL T X 



X = 
Y = 



Y 



"tt~ t A' 

U l A \2 





'A' 

A \2 




A' 
.22. 




'A' 

A \2 




A' 

.22. 


A 1 
A \2 




A' 

^22 





A'22 U 2 U 1 A\ 



12 
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The computation of B x is straightforward from the equation TB 1 = B 

B 1 = L- T U^ T PB. 
We are now able to write the algorithm. 



Algorithm 5.1 Kalman form 



Require: Aan x n matrix over a field, B, a n x m matrix 
Ensure: r, T, H, X, Y, B x as in theorem 4.1 

1 : (V, r) = CompressedKrylovMatrix(j4, B) 

if (r=n) then 

return (n, Id, A, 0, 0, B) 
else 

(L,[C/ 1 C/ 2 ],P) = LUP(0 



9 
10 
11 
12 
13 
14 
15 
16 



T = 

B 1 

A' - 
X -- 



V 



P 1 







L- T U^ T PB 



PAP T 



A' A 1 

A' 



L- T U^ T A' 12 
A' 



21 



12 

1' 

i 22 



^ ~ ^22 ^2 ^1 ^12 

for all j do 

m,j = IjLJ 1 as explained in section 5.1 
end for 

Build the polycyclic matrix H using the mj as shown in section 5.1 . 
return (r t T,H,X,Y,B{) 
end if 



Lastly, note that the LUP factorization of V T is already computed at the end of the call to 
CompressedKryiovMatrix. Thus, step 5 in algorithm 5.1 can be skipped. 

5.3 LU-Krylov : a cubic variant 

In [DPW05], we introduce an algorithm for the computation of the characteristic polynomial : LUK. 
Alike Keller-Gehrig's algorithm, it is also based on the Krylov iterates of several vectors and relies 
as much as possible on matrix multiplication. But the krylov iterates are computed with matrix 
vector products, so as to avoid the logn factor in the time complexity. As a consequence it is 
0(n 3 ) algorithm, but we showed that it was faster in practice. 

Algorithm 5.2 shows how to adapt this algorithm to the computation of the Kalman form. We 
expect this algorithm to be the more efficient in practice. 
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Algorithm 5.2 Kalman-LUK : Kalman-LU-Krylov 



Require: Aan x n matrix over a field, B, a n x m matrix 
Ensure: r, T, i7, X, Y, S x as in theorem 4.2 

1 : u = Bi 

K == Av A^v 

'") T \ „ {The matrix A is computed on the fly : at most 



(L, [Ui\U 2 ],P) = LUP(A T ),n = rank(A) 
2ri columns are computed} 

3: m = (mi, . . . , m r i) = L^+iL^ 1 ^ 
4: / = - Elii miX*" 1 

5: if (n = n) then 
6: return (n, Id, A, 0, 0, _B) 
7: else 

8: A' = PAP T 



9: 
10: 

11: 

12: 
13: 
14: 

15: 

16 
17 
18 

19: 

20: 

21: 

22: 

23 
24 
25 



A 1 A' 
A 1 A' 

A 21 ^22 



where A' n is n x n. 



B' = 



L- T 
/ 



-u?ur T 



PB 



Compute the permutation Q s.t. B'Q 
if (/i = 0) then 



{Z is (n - n) x /i} 



y : 
T = 



A 







return {n,T,C f ,X,Y,X) 
else 

(r 2) T (2) , # (2) , X<- 2 \ YV\ B{ 2 } ) = Kalman-LUK(A fl , Z) 




T = 



A 



T( 2 ) 



J = L- T tfr T Ai 2 T( 2 ) = [Ji J 2 ] {Ji is n x r 2 and J 2 , n x (n - n - r 2 )} 
" " o i/( 2 ) 

J 2 

return (n +r 2 ,T, Y( 2 ),Bi) 
end if 
end if 
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